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<"^ ' Abstract 

Harmonic polylogarithms H(a;a;), a generalization of Nielsen's polylogarithms S„, p (x), appear fre- 
<**> |_ quently in analytic calculations of radiative corrections in quantum field theory. We present an algorithm 

for the numerical evaluation of harmonic polylogarithms of arbitrary real argument. This algorithm is 
<-^h implemented into a FORTRAN subroutine hplog to compute harmonic polylogarithms up to weight 4. 
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PROGRAM SUMMARY 

Title of program: hplog 

Version: 1.0 

Release: 1 

Catalogue number: 

Program obtained from: Thomas .GehrmannScern. ch, Ettore.Remiddi@bo.infn.it 

E-mail: Thomas. Gehrmann@cern.ch, Ettore.Remiddi@bo.infn.it 

Licensing provisions: none 

Computers: all 

Operating system: all 

Program language: F0RTRAN77 

Memory required to execute: Size: 516k 

No. of lines in distributed program: 3271 

Other programs called: none 

External files needed: none 

Keywords: Harmonic polylogarithms, Fcynman integrals 

Nature of the physical problem: Numerical evaluation of the harmonic polylogarithms up to weight 4 for 
arbitrary real argument. These functions are emerging in Fcynman graph integrals involving more than 
one mass scale. 

Method of solution: For small values of the argument: scries representation; other values of the argument: 
transformation formulae. 

Restrictions on complexity of the problem: limited up to HPL of weight 4, the algorithms used here can be 
extended to higher weights without modification. 

Typical running time: On average 0.3 ms for the evaluation of all harmonic polylogarithms up to weight 4 
on a Pentium III/600 MHz Linux PC. 



LONG WRITE-UP 



1 Introduction 

The Euler Dilogarithm 1^2(2;), and its generalizations, Nielsen's poly logarithms O, have been playing a 
central role in the analytic evaluation of integrals arising in perturbative quantum field theory. A reliable 
and widely used numerical representation of these functions (GPLOG) [|2| has been available for already thirty 
years. Going to higher orders in perturbation theory, it was realized recently that Nielsen's polylogarithms 
are insufficient to evaluate all integrals appearing in Feynman graphs at two loops and beyond. This 
limitation can only be overcome by the introduction of further generalizations of Nielsen's polylogarithms, 
the Harmonic Polylogarithms (HPLs). The HPLs, introduced in j|, together with their "2-dimensional" 
extensions m , are already now playing a central role in the analytic evaluation of Feynman graph integrals H, 
pj. HPLs appear also as inverse Mellin transformations of harmonic sums, which were investigated and 
implemented numerically in H. 

Let us recall here two key features of an analytical calculation. The first is to express everything in 
terms of mathematical quantities whose properties are known (so that in particular the final formulae can 
be written in a unique way, which does not hide cancellations between the various terms involved). The 
second is to allow the fast and precise numerical evaluation of all the mathematical quantities introduced. 
Obviously, such fast and precise evaluation relies on the knowledge of the analytical properties. 

In this paper, somewhat as a continuation of dealing with the numerical evaluation of Nielsen's 
polylogarithms M (a subset of the HPLs), we shortly review the analytical properties of the HPLs and 
then show how they can be used for writing a FORTRAN code which evaluates the HPLs up to weight 4 (see 
Section below for the definition of the weight; 4 is the maximuum weight required in the calculations 
of W) with absolute precision better than 3 x 10~ 15 (i.e. standard double precision) with a few dozens 
of elementary arithmetic operations per function (just a single dozen in the most favourable case). Given 
the large number (81) of HPLs of weight 4 and the many algebraic relations among them, our FORTRAN 
routine evaluates the whole set of all the HPLs up to the required weight - at variance with GPLOG M, 
which evaluates separately the various Nielsen's polylogarithms. 

The plan of the paper is as follows. Section recalls the definitions of the HPLs. Their algebraic proper- 
ties are discussed in Section 0, where we show how to use these properties for separating the functions into 
reducible and irreducible ones. Section [| studies the analytic properties which allow to perform converging 
power series expansions and to accelerate their convergence. Relations between HPLs for different ranges 
of the arguments are derived in Section |5|. Section |6| explains how the properties recalled above are used to 
implement the HPLs into a FORTRAN subroutine hplog, and Section |7| how the correct implementation is 
checked. Finally, we describe the usage of the subroutine hplog in Section |3| and provide a few numerical 
examples in Section |9|. 

2 Definitions 

The Harmonic Polylogarithms (HPLs), introduced in ||], are one- variable functions H(a;x) depending, 
besides the argument x, on a set of indices, grouped for convenience into the vector a, whose components 
can take one of the three values (1, 0, —1) and whose number is the weight w of the HPL. More explicitly, 
the three HPLs with w = 1 are defined as 

f x dx' 
H(l;aO = / -^- 7 = -ln(l-x), 
Jo l-# 

H(0; x) = \nx , 
H(-l;z) = f X J^-= H l + x); (2.1) 

In 1 -h X 



their derivatives can be written as 



dx 
where the 3 rational fractions f{a;x) are given by 



— E(a;x)=f(a;x) , a = 1,0,-1 (2.2) 



f(l;x) = 
f(0;z) = 



1 



1-x ' 
1 



X 

f(-l;a:) = -i-. (2.3) 

1 + x 

For weight w larger than 1, write a = (a, b), where a is the leftmost component of a and b stands for the 
vector of the remaining (w — 1) components. The harmonic polylogarithms of weight w are then defined 
as follows: if all the w components of a take the value 0, a is said to take the value 0^, and 

R(6 w ;x) = ^rhx w x , (2.4) 



while, if a ^ (\ 



H(a;x) = / da;' i(a;x') E(b;x') . (2.5) 

Jo 



In any case the derivatives can be written in the compact form 

— H(a; x) = f(a; x)B(b; x) , (2.6) 

da; 

where, again, a is the leftmost component of a and b stands for the remaining (w — 1) components. 



It is immediate to see, from the very definition Eq. (2.5), that there are 3 10 HPLs of weight w, and that 
they are linearly independent. The HPLs are generalizations of Nielsen's polylogarithms B. The function 
Sn tP (x), in Nielsen's notation, is equal to the HPL whose first n indices are all equal to and the remaining 
p indices all equal to 1: 

S ntP (x) = B.(6 n X;x); (2.7) 

in particular the Euler polylogarithms Li„(a;) = Sn-i^x) correspond to 

Li„(x)=H(0 n _i,l;a;). (2.8) 

All relations between HPLs which will be used for their numerical evaluation are easily derived by using 



recursively the definition (2J3) and the related differentiation formula (2.6) 



3 The algebra and the reduction equations 

As shown in H , the product of two HPLs of a same argument x and weights p, q can be expressed as a 
combination of HPLs of that argument and weight r = p + q, according to the product identity 

R(p;xMq;x) = £ K(r;x) , (3.1) 

r—p^±iq 

where p, q stand for the p and q components of the indices of the two HPLs, while pttl q represents all 
mergers of p and q into the vector r with r components, in which the relative orders of the elements of p 
and q are preserved. 

The explicit formulae relevant up to weight 4 are 

H(a; x) H(6; x) = H(a, 6; x) + H(&, a; x) , (3.2) 



H(a;a;) H(&, c; x) = H(a, 6, c; a;) + H(fo, a, c; x) + H(&, c, a; a;) , 
H(a;x) H(&, c, d;x) = H(a, 6, c, d; x) + H(6, a, c, d; x) + H(6, c, a, d; x) + H(6, c, d, a; a;) , (3.3) 

and 

H(a, 6; a;) H(c, d; x) — H(a, 6, c, d; x) + H(a, c, b, d; x) + H(a, c, d, 6;x) 

+ H(c, a, 6, d; x) + H(c, a, d, 6; a;) + H(c, d, a, 6; x) , (3-4) 

where a, 6, c, d are indices taking any of the values (1, 0, —1). The formulae can be easily verified, one at 
a time, by observing that they are true at some specific point (such as x — 0, where all the HPLs vanish 
except in the otherwise trivial case in which all the indices are equal to 0), then taking the a;-derivatives of 
the two sides according to Eq. (|2.6[ ) and checking that they are equal (using when needed the previously 
established lower- weight formulae). 



Another class of identities is obtained by integrating (2.4) by parts. These integration-by-parts (IBP) 
identities read: 

if (mi, . . . , m q ; x) = H{m\; x)H(m,2, . ■ ■ , m q ; x) — H(rri2, mi; x)H(rri3, . . . , m q ; x) 

+ ... + {-lf +1 H{m q ,...,m 1 -x). (3.5) 

These identities are not fully linearly independent from the product identities. 



By using Eqs.(3.2) for all possible independent values of the indices a, b one obtains 6 independent 
relations between the 9 HPLs of weight 2 and the products of 2 HPLs of weight 1 ; those relations can be 
used for expressing 6 of the HPLs of weight 2 in terms of 3 HPLs of weight 2 and products of 2 HPLs of 
weight 1. The choice of the 3 HPLs (referred to, in this context, as "irreducible") is by no means unique; 
by choosing as irreducible HPLs of weight 2 the 3 functions H(0, 1; a;),H(0, — 1; x) and H(— 1, l;x), the 
reduction equations expressing the 6 "reducible" HPLs of weight 2 in terms of the irreducible HPLs read 

H(l,l;x) = hi(l;x)R(l;x) , 

H(l,0;x) = H(l;a;)H(0;x)-H(0,l;a;) , 
H(l,-l;x) = H(l;ir) H(-l; a;) - H(-l, 1; x) , 

H(0,0;») = -H(0;»)H(0;a:) , 

H(-l,0;x) = H(-l;a;)H(0;»)-H(0,-l;a;) , 

H(-1,-1;jc) = -H(-l;*)H(-l;x) . (3.6) 

Similarly, at weight 3 one has 27 HPLs and 19 independent product and integration-by-parts identities, 
expressing 19 reducible HPLs in terms of 8 irreducible ones; at weight 4 there are 81 HPLs, 63 independent 
identities, and correspondingly 63 reducible and 18 irreducible HPLs. 

4 The analyticity properties 

At weight 1, the HPLs are just logarithms; with the standard conventions for the specification of the 
logarithm, it is immediately seen that H(l; x) has a cut along the real axis from x = 1 to x — > +oo, H(0; x) 
from x = to x — ► — oo and H(— l;i) from x — —1 to x —> — 00. With the usual +ie prescription the 
corresponding complex values along the whole real axis are given by 

R(l;x + ie) = -ln(|l-a:|)+wr0(a;-l) , 
R(0;x + ie) = ln(|a;|) + wr0(-a:) , 
H(-l;ar + ie) = lnfll + x\) +iir6(-x - 1) . (4.1) 

In particular H(l; x) and H(— 1; x) are analytic at x — and can therefore be expanded in series of powers 
of x around x — 0, the first term of the expansion being of order x. 



At weight 2, let us start from the HPLs whose rightmost (trailing) index is 1, 

R(a,l;x) = dti(a,t)R{l:t) , a =1,0,-1. 

Jo 

For a = 1, i?(l,l;x) has the same cut from x — 1 to x — +oo as H(l;x). For a — 0, by recalling that 
H(l;£) can be expanded in powers of t at t = 0, and that the first term of the expansion is t {i.e. the 
constant term vanishes), one finds that also 



H(0,l;aO= / ^m;t) 
Jo l 



shares the same analyticity properties as H(l; x), namely it has the same cut from x = 1 to x = +oo and 
it can also be expanded in powers of x around x = 0, the first power of the expansion being x (i.e. the 
constant term vanishes). For a = — 1, finally, one sees that 



K(-l,l;x) = l X ^- t K m 



can again be expanded in x at x = (the first power of the expansion being in this case a; 2 ); but besides 
the right cut 1 < x < +oo implied by the presence of the H(l;t) in the definition, H(— 1, l;x) has also the 
left cut — oo < x < — 1 due to the 1/(1 + 1) fraction. The two cuts can be easily separated by writing 



H(-l,l;a;) = H+(-1,1;.t)+H_(-1,1;x) , 

r dt 

H + (-l,l;a;) = jf — pH(l;*)-H(l;-l)] , 



H-(-l,l;aO = £^^(15-1)] , (4.2) 

where H + (— 1, 1; x) has only the right cut 1 < x < +oo and H_(— 1, 1; x) only the left cut — oo < x < — 1, 
(and both admit an expansion in x at x = whose constant term vanishes). Note that in the above example 
the value of H(l; —1) = — ln2 is finite, as H(l; x) is regular at x = — 1. 

The same discussion applies to the functions H(a, — l;x),a = 1,0,-1; H(— 1,— l;x) and H(0, — l;x) 
share the same analyticity properties as H(— 1; x), while 



H(l,-l;x) = £ T ^H(-l 



-1-t) 

develops, besides the left cut — oo < x < — 1 implied by H(— l;t), the right cut 1 < x < +oo due to the 
fraction 1/(1 — t) in the definition, for which one can write: 



H(l,-l;a;) = H+(l,-l;a;)+H_(l,-l;x), 
H+(l,-l;x) = J JL-K(-l-l), 

f' x A+ 

H_(l,-l;x) = / — — [H(-l;t)-H(-l;l)] , (4.3) 

Jo i — t 

where H(— 1; 1) = ln2 is well defined as H(— 1; x) is regular at x — — 1. 

The procedure might be extended to cover also the HPLs whose rightmost index is 0. It turns out, 
however, that it is not really necessary to work out explicitly also this last case (which is more complicated, 
as also the — oo < x < cut is present). One can in fact exploit the reduction formulae of the previous 
section for choosing a set of irreducible HPLs with the rightmost index always different from 0, (with the 
exception of the w = 1 case, where the irreducible set coincides with the set of all 3 HPLs), and then 
express all the HPLs with w > 1 and rightmost index in terms of those irreducible HPLs. 



It is clear from the above discussion that any HPL of arbitrary weight, whose rightmost index is I or 
— I and whose other indices are otherwise arbitrary, say H(a;x), can always be separated into the sum of 
two functions, 

H(a;x) =H+(a;a;)+H_(a;.T) , (4.4) 

where H+(a;x) has only the right cut I < x < +00 and H_(a;x) only the left cut —00 < x < — 1. (The 
formula applies of course even when one of the two cuts is missing, as in that case the corresponding 
function is equal to zero; one can write for instance H(0, 1, 1; x) = H + (0, 1, 1; a;), as H_(0, 1, 1; x) — 0.) All 
the involved functions (the HPLs as well as the functions corresponding to the separated cuts with branch 
points at 1 or — 1) are regular at x = and can be expanded in series of powers of x, the constant term 
always being zero. 

Algorithmically, if a — (a,b) (i.e. a is the leftmost index of a, and b stands for the remaining w — 1 
indices) and H(&; x) admits the separation of the cuts 



H(6;x) =H+(6;x)+H_(6;a;) 



(4.5) 



with the two functions H + (6; x) and H_(&; x) corresponding to the right 1 < x < +00 and left —00 < x < —I 
cuts, the explicit separation formulae for H(a, 6; x) of weight w + 1 



H(a, 6; x) 
in the three cases a = 1, 0, — 1 are given by 



R+(l,b;x 



H_(l,6;a; 

H+(0,&;x 

H_(0,6;x 

H+(-l,fe;.T 

H_(-l,6;x 



H + (a, 6; x) + H_(a, b; x) 



H+(6;a;)+H_(6;l) 
H_(?;a;)-H_(6;l) 



dt \ 


1-t i 
dt \ 


1-t 1 



x dt 



H+(6;x) 
-H_(6;x) 



dt 

„ l + t 
dt 



R+(b;x)-R + (b;-l) 

H_(6;x)+H + (6;-l) 



/o l + < 
The separation of the two cuts is illustrated with two explicit examples: 

H(-l,l,l;y) = H + (-l,l,l;y) + H_(-l,l,l;y) 
with: H+(-l,l,l;y) = H(-l, 1, 1; y) - I ln 2 2 H(-l; y) ., 

H_(-l,l,l;y) = iln 2 2H(-l;y), 

and 



(4.6) 



(4.7) 



H(0,l,-l,l;s) = H+(0, 1,-1, l;a;)+H_(0, l,-l,l;a:), 

with: H+(0, 1,-1,1; x) = H(0, 1, -1, l;x) - ln2 H(0, 1, -l;x) + ln 2 2 H(0, l;x) , 

H_(0, 1,-1,1; x) = ln2H(0,l,-l;a;)-ln 2 2H(0,l;x) . (4.8) 

The actually chosen irreducible HPLs are listed in Table [y. Note that, in addition to the requirement 
of the absence of the index in the rightmost place, one can also imp ose the absence of the index 1 in 
the leftmost place. That guarantees, owing to the definition Eq. (|2.5| ), that all the HPLs of the above 
irreducibe set have finite values at x = 1. The explicit values at x — 1 are needed in the separation of cuts, 



as well as in some of the transformation formulae discussed in the following section (as a matter of fact, 
the x = 1 values can also be obtained as consistency conditions of the many transformation formulae of 
various kinds, which can be easily established for the HPLs). Up to weight 4, they can easily be obtained 
by using the table of definite integrals given in |7J , and we list them in Table [y as well. It can be seen 
that, besides ln.2 = 0.693147180559945 . . . and tt 2 = 9.869604401089359 . . ., only two more transcendental 
constants appear: C 3 = 1.202056903159594 ... and Li 4 (l/2) = 0.517479061673899 . . .. The values at x = -1 



can be obtained directly from the transformation discussed in Section 5.1 below 



5 Identities for HPLs of related arguments 

We discuss here the identities valid for HPLs whose arguments are related by the transformations x = — y, 
x = 1/t and x = (1 - r)/(l + r). 

5.1 The x = —y transformation 



For real x (in the whole range — oo < x < +oo), if y = —x, at weight w = 1 one has from Eq. (4.1) 



H(l;s + *e) = -H*(-l;i/ + »e) , 
H(0;x + ie) = ( R(0,y + ie) -m )* , 
H(-l;a; + *e) = -H*(l;y + te) , (5.1) 

where the asterisk stands for the complex conjugate. 

At w > 1, if the rightmost index is different from (as is the case for the irreducible HPLs), writing 
explicitly the indices up to w — 4 one finds immediately the following formulae, which apply for a\ = 1,-1 
and any value of the other indices: 

E{a 2 ,av,x + ie) = (-l) ai+Q2 H*(-a 2 , -any + ie) , 
K(as,a 2 ,av,x + ie) = (-l) ai+a2+a3 H*(-a 3 , -a 2 , -ar,y + ie) , 
K(a 4 ,a 3 ,a 2 ,av,x + ie) = (-l) ai+a2+a3+a4 H*(-a 4 , -a 3 , -a 2 , -any + ie) . (5.2) 

5.2 The x = 1/t transformation 

For real x in the range 1 < x < +oo, put 



1 1 

t x 



l<x<+oo, 1 > t > . (5.3) 



At weight 1, one has 



H(l;a; + ie) = H(l;t) +H(0,t) +in , 
R(0;x) = -H(0,t) , 
H(-l;a;) = H(-l;t) -H(0,t) . (5.4) 

For higher weight, one can proceed recursively (i.e. by induction on the weight w), observing further that 
for w > 1 it is sufficient to establish the identities for the irreducible HPLs only, as the other cases can be 
obtained through the reduction formulae. Let us consider a generic irreducible HPL, say H(a, b\x = 1/t), 
and define accordingly 

X(a,b;t) =H(o,6;l/*) . 

Quite in general, 

X(a, b; t) = X(a, b;l)+ / dt'—X(a, b; t') ; 
./i at 



by using, in the r.h.s., the very definition of X(a,b;t) and Eq. (2.5) the previous equation reads 

H L, 6; i) = H(o, b; 1) - J dt'± f (a, £) H (V; ~) . 

It was already remarked that the irreducible HPLs can be chosen with the leftmost index different from 1 ; 
therefore in the above equation the index a can take only the values and — 1; correspondingly one finds 



h(o,&;M = R(0,b-,1)-J dt'-H^i 



H(-M;i) = H(-U;l)-| d t> (i-^) H(^I) . (5.5) 

By proceeding recursively from lower to higher weights, H(b; 1/t') can be considered as already given in 
terms of HPLs of argument t', so that due to Eq. (|]q) X(a,b;t), hence K(a,b;x = 1/t) is also expressed 
in terms of HPLs of argument t. 

As an example of this transformation, one finds for instance, 

H(-l, 1; x) = -H(0, 1; t) - H(0, -1; t) + H(-l, 1; t) - ^H(0; t)H(0; t) + H(-l; i)H(0; t) + ^ 



+in{ H(-l;i)-H(0;i) - ln2 

H(0, 1, 1; x) = H(0, 0, 1; t) - H(0, 1, 1; t) + H(0; t) f y - H(0, 1; t) - ^H(0; t)H(0; t)) + Cs 

+*tt f-H(0, 1; t) - ^H(0; i)H(0; *) + y ) ■ (5-6) 

5.3 The x = (1 — r)/(l + r) transformation 

The transformation 

x = , r = , < x < +oo , 1 > r > -1 . (5.7) 

1 + r 1 + x 

can be treated similarly. For x in the range < x < 1, r is positive, 1 > r > 0, and all the involved 
functions are real, while for x in the range 1 < x < +oo the variable r is negative, > r > — 1; at weight 
w = 1, therefore, one has 

H(l;a:) = -H(O;|r|)+H(-l;r)-ln2 + *7r0(-r) , 
H(0;z) = -H(l;r) -H(-l;r) , 
H(-l;x) = -H(-l;r)+ln2. (5.8) 

To work out the corresponding formulae for any of the irreducible HPL of higher weight, say H(a, b;x = 
(l-r)/(l + r)), define 

Y{a,b-r)=n(a,b-\-^ 
V 1+r 

and evaluate it through the obvious relation 

Y(a,b:r) = Y{a,b;0) + / dr'-— r(o, b;r') 



o dr' 



which on account of the definitions becomes 



II ( a, b; L-l) = H(o, 6; 1) - / dr' — ^— -f fa, I— ^ H (b; * ' ' 



1 + r J v ' ' ' 7 (1 + r') 2 V ' 1 + ? - '/ V ' ! + r ' 



It was already recalled that the leftmost index of an irreducible HPL does not take the value 1; for a — 0, — 1 
an elementary explicit calculation gives 

Proceeding again recursively from lower to higher weights, H(6; (1— r')/(l+r')) can be considered as already 



given in terms of HPLs of argument r'; therefore due to Eq. (2.5) Y(a,b;r), hence H(a,6;x = (1 — r)/(l+r)), 
is also expressed in terms of HPLs of argument r. 

As an example of the above identities, written in a form that applies in the whole range < x < +oo 
one finds for instance 

H(0, 1; x) -H(0, 1; r) - H(0, -1; r) + H(-l, 1; r) 

H(l;r) + H(-l; r) ) ( H(0; |r|) + ln2 - wr0(-r 



_^H(l;r) + iH(-l;r))H(-l;r) + l-, 

H(-l,l;x) = -H(O,-l;r)+H(-l;r)(H(O;|r|)+m2-wr0(-r) 

I , s , s 7T 2 ln 2 2 

--H(-l; r)H(-l; r) + — - — . (5.10) 



As a last remark concerning the transformation (5.7), let us observe that its fixed points satisfy the 

equation 

1 - x 

x= Tl~ ' 
1 + x 

whose solutions are 

x+ = + (v/2 - 1) , 

x_ = -(V2 + 1). (5.11) 

6 The numerical evaluation 

The analytical properties discussed in the previous sections can be exploited to get a fast and precise 
numerical evaluation (absolute error less than 3 x 10~ 15 ) of the HPLs of real but otherwise arbitrary 
argument x. We divide the whole range — oo < x < +oo into the 5 regions 

• -oo < x < -(\/2+ 1), 

• -(n/2 + 1) <x<-(v^-1), 
. -(\/2-l) <a:<+(V2-l), 

• (\/2-l) <x< (\/2 + l), 

• (\/2 + 1) <x< +oo. 

In evaluating all HPLs, it is always assumed that the (infinitesimal) imaginary part of the argument is 
positive, i.e. the HPLs are evaluated for argument x + ie. 

All transformation and expansion formulae used in the numerical evaluation of the HPLs were obtained 
using the computer algebra program FORM S . 



6.1 The central region -{y/2 - 1) < x < +( v / 2 - 1) 

In the central region — (\/2 — 1) < x < +(v / 2 — 1), we evaluate H(0; x), according to its definition Eq. (pH|), 
as ln(|a;|) + iTr9(—x), and all the other irreducible HPLs up to the required weight through a power series 
expansion in x around x — 0. 

As |a;| < (-\/2 — 1) = 0.4142 . . ., to obtain the aimed 3 x 10~ 15 absolute precision, one should keep in 
principle 34 terms in the expansion in powers of x. If the indices of some HPLs take only the values and 
1, it has already been observed that the corresponding function has only the right cut 1 < x < +oo; in that 
case, the convergence of the power series expansions improves considerably by the Bernoulli H change of 
variable 

x = l-e~ u , u = - ln(l - x) . (6.1) 

This transformation is invertible in the strip \^su\ < 2ir and moves the x = 1 branch point to u = +oo. By 
expanding in powers of u at u = one obtains a series in u whose radius of convergence is 2ir (as opposed 
to 1 for the series in x). As x moves from — (V2 — 1) through to {s/2 — 1), u/(2tt) moves from —0.055 . . . 
through to 0.085..., so that in the worst case |u/(27r)| ~ 0.085...; to attain a 10 -16 precision it is 
sufficient to keep up to 15 terms in u. A few terms drop when using standard Chebyshev economization, so 
that in practice the required accuracy is obtained with at most 12 Chebyshev polynomials in u. The same 
approach works when the indices of the HPLs take only the values and — 1 ; in that case the corresponding 
function has only the left cut — oo < x < — 1 and the appropriate change of variable is 
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v = ln(l + x) 



(6.2) 



it is found again that after Chebyshev economization at most 12 Chebyshev polynomials in v arc needed. 
When 1 and —1 both appear as indices of the HPLs both cuts are present. An expansion in powers 
of x at x — can be equally carried out, but the above changes of variable do not help for speeding up 
convergence. Indeed, by using for instance the transformation Eq. (pTj), the left branch point at x = — 1 
is mapped into v — — ln2, so that the radius of convergence of the u expansion at u = would be in fact 
reduced from 1 to ln2 ( ln2 = 0.6931 . . .). This problem is, however, easily overcome by using Eq. (4.4) and 



evaluating independently, for each irreducible H(a;a;), the two related functions H + (a;x) and H_(a;a;) by 
using the variable u of Eq. (6.1) for H+ (a; x) and the variable v of Eq. (6.2) for H_(a; ir). It should be noted 
that it is not actually necessary to evaluate both H + (a; x) and H_(a; x) for each irreducible H(a; x). Indeed, 
if the rightmost index is 1, H + (a;x) differs from H(a ; x) by a com bination of constants times HPLs with 
the left cut and smaller weight, see for instance Eqs. (4/7) and ( [4.8| ), while H_(a;cc) is just that difference. 
In a systematic approach to the numerical evaluation of the HPLs in order of increasing weight, such lower 
weight functions can be considered as already evaluated. Therefore only H + (a; a;) has to be evaluated from 
scratch as a power series in the suitable variable. The obvious analogue applies when the rightmost index 
is -1. 

Some explicit examples of actual expansions in u and v follow: 
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The fast convergence of the expansion is shown by the fast decrease of the higher order coefficients. In the 
first and third of the above equations the fast decrease is manifest; in the second equation it comes from the 
strong cancellations between the two terms in brackets; for instance, the actual value of the coefficient of 
the w 4 term is 9.09 . . . x 10~ 4 , the first rational fraction being 3/4 (the coefficient could indeed be used for 
obtaining the rational approximation 9/13 = 0.692307... to ln2 = 0.6923147...), and the effect is of course 
enhanced in the subsequent terms. Some care has therefore to be taken in implementing these coefficients 
in a numerical program. To avoid large cancellations inside the coefficients during the execution of the 
program, we converted all coefficients from combinations of rational fractions and transcendental constants 
into real numbers of the desired accuracy, so that in the FORTRAN code each power of u and v is multiplied 
by a single constant in double precision. 

It should be noted that u and v are both in the range [— ln(v2 + 2) : ln(v2 + 2)]. The Chcbyshev 
economization is carried out on those terms by rescaling u and v with a factor 20/11, to an interval slightly 
smaller than [—1 : 1] . Using this only approximate rescaling allows an extension of the numerical formulae 
used in this region also slightly beyond its boundaries, which will be used as a check on the numerical 
program below. 

6.2 The region (V2-1) <x < (V2 + 1) 

In the region (\/2— 1) < x < (V2+ 1) we use the change of variable of Eq. (5/7) in Section 5JJ for expressing 



the irreducible HPLs of argument x in terms of the irreducible HPLs of argument r, with 



1 



Note that x = (v2 — 1) corresponds to r — (v2 — 1) (one of the fixed points of the transformation), while 
x = (v2+ 1) corresponds to r = — (v2~ — 1), so that the region (v2 — 1) < x < (v2+ 1) is mapped into the 
r-range — (v2— 1) < T < (v2 — 1), which is exactly the central region discussed in the previous subsection. 
We evaluate the HPLs of r with the series expansion discussed above, then obtain the required HPLs of x 



by means of the transformation formulae discussed in Section 5.2 
6.3 The region (\/2 + 1) < x < +oo 



We use the change of variable x = 1/t of Eq. (5.2), discussed in Section 5.2, for expressing the irreducible 
HPLs of argument x in terms of the irreducible HPLs of argument t, with t — 1/x. Note that x = (\/2~ + 1) 
corresponds to t = (v2— 1), while x — > +oo translates into t — ► 0, so that the region (v2+ 1) < X < +oo is 
mapped into the t-range < t < (\/2 — 1), which is the positive half of the central region discussed at the 
beginning of this section. We evaluate the HPLs of t from their series expansion, then obtain the required 



HPLs of x by means of the transformation formulae of Section 5.2 



6.4 The region -(-v/2 + 1) <x< -(V2-1) 



We use the change of variable x = —y. We evaluate the HPLs of argument y = —x as in Subsection 6.2 



and then use the formulae of Subsection 5.1 for expressing the required HPLs of argument x in terms of 
the HPLs of argument y. 

6.5 The region — oo < x < — (\/2 + 1) 



Essentially the same as before. We evaluate the HPLs of argument y = —x as in Subsection |6.3| and then 



use the formulae of Section 5.1 
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7 Checks 

We have carried out several checks on our implementation of the algorithm described in the previous section 
into a FORTRAN subroutine. 

An immediate check of the numerical implementation of the HPLs is provided by the derivative for- 



mula Eq. (2.6). Evaluating the left-hand side of Eq. fl2.q ) numerically with a standard symmetric 4-point 
differentiation formula, and comparing it with the right-hand side evaluated directly, we found agreement 
within an accuracy of 10 -12 or better. This accuracy is mainly limited by rounding errors induced by 
the small interval size used in the differentiation formula. We used 10 as interval size, which implies a 
theoretical accuracy of about 10 -16 . This accuracy is however reduced to the observed 10 -12 by rounding 
errors arising from taking the difference between the function values evaluated at interval points spaced by 
10~ 4 only. 

We also checked the continuity of all HPLs across the boundaries of the different regions introduced 
in Section 0, which match onto each other at the points ±x + and ±ai_ of Eq. (pMl). Evaluating the 



HPLs according to the algorithms appropriate to the regions left and right of the boundaries, we found 
both limiting values to agree within 3 x 10~ 15 or better for both real and imaginary parts. Moreover, 
we evaluated the HPLs in a few points scattered in a small neighbourhood (of size ±10 -10 ) across the 
boundaries of the regions discussed in the previous section, by using for each point the two different 
algorithms used separately in the FORTRAN code at each side of the boundaries and then comparing the 
results. Again, we found agreement within the desired accuracy of 3 x 10 5 . 

8 The subroutine hplog 

8.1 Syntax 

The routine hplog has the following syntax: 

subroutine hplog (x,nw, He 1 ,Hc2,Hc3,Hc4, 
$ Hrl,Hr2,Hr3,Hr4,Hil,Hi2,Hi3,Hi4,nl,n2) 

8.2 Usage 

In calling hplog, the user has to supply 

x: The argument for which the HPLs are to be evaluated, x is of type real*8. It can take any real 
value — oo < x < oo. 

nw: The maximum weight of the HPLs to be evaluated, nw is of type integer. It is limited to 1 < nw < 4. 

nl ,n2: Define the indices for which the HPLs are evaluated, (nl ,n2) are of type integer. Allowed combina- 
tions are (—1, 1) (evaluate all HPLs), (0, 1) (evaluate only HPLs with and 1 appearing in the vector 
of arguments) and (—1,0) (evaluate only HPLs with and —1 appearing in the vector of arguments). 

The output of hplog is provided in the arrays He 1 , Hc2 , Hc3 , Hc4 , Hr 1 , Hr2 , Hr3 , Hr4 , Hi 1 , Hi2 , Hi3 , Hi4. 
These have to be declared and dimensioned by the user as follows: 

complex*16 Hcl ,Hc2,Hc3,Hc4 

real*8 Hrl ,Hr2,Hr3,Hr4 

real*8 Hil ,Hi2,Hi3,Hi4 

dimension Hcl(nl:n2) ,Hc2(nl :n2,nl :n2) ,Hc3(nl :n2,nl :n2,nl :n2) , 
$ Hc4(nl:n2,nl:n2,nl:n2,nl:n2) 

dimension Hrl(nl:n2) ,Hr2(nl :n2,nl :n2) ,Hr3(nl :n2,nl :n2,nl :n2) , 
$ Hr4(nl:n2,nl:n2,nl:n2,nl:n2) 

dimension Hil(nl:n2) ,Hi2(nl :n2,nl :n2) ,Hi3(nl :n2,nl :n2,nl :n2) , 
$ Hi4(nl:n2,nl:n2,nl:n2,nl:n2) 
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It should be noted that this declaration is always needed, even if hplog is called with nw< 4. After calling 
hplog for a given argument x, the complex arrays Hcl, Hc2, . . . contain the complex values of the 
corresponding HPLs of weight 1,2,..., the real arrays Hrl , Hr2 , ... contain the real values and Hil , 
Hi2, . . . the immaginary parts divided by n, so that, for instance, the following equality holds 

Hc3(kl,k2,k3) = cmplx( Hr3(kl ,k2,k3) , pi*Hi3(kl ,k2,k3) ) 

The complex parts of the HPLs are always evaluated assuming an infinitesimal positive imaginary part of 
the argument. 

The subroutine does not need initialization. 

8.3 Example 

The following example program illustrates how to evaluate HPLs up to weight 2 for a given value of x, and 
to write out their real parts: 

program testhpl 

integer nl ,n2,nw, il , i2 

parameter (nl=-l) 

parameter (n2=l) 

parameter (nw=2) 

complex*16 Hcl ,Hc2,Hc3,Hc4 

real*8 Hrl ,Hr2,Hr3,Hr4 

real*8 Hil ,Hi2,Hi3,Hi4 

real*8 x 

dimension Hcl(nl:n2) ,Hc2(nl :n2,nl :n2) ,Hc3(nl :n2,nl :n2,nl :n2) , 
$ Hc4(nl:n2,nl:n2,nl:n2,nl:n2) 

dimension Hrl(nl:n2) ,Hr2(nl :n2,nl :n2) ,Hr3(nl :n2,nl :n2,nl :n2) , 
$ Hr4(nl:n2,nl:n2,nl:n2,nl:n2) 

dimension Hil(nl:n2) ,Hi2(nl :n2,nl :n2) ,Hi3(nl :n2,nl :n2,nl :n2) , 
$ Hi4(nl:n2,nl:n2,nl:n2,nl:n2) 

write(6,*) 'Input x: ' 

read(5,*) x 

call hplog (x,nw, He 1,Hc2,Hc3,Hc4, 
$ Hrl,Hr2,Hr3,Hr4,Hil,Hi2,Hi3,Hi4,nl,n2) 

do il = nl,n2 

write(6,101) il,Hrl(il) 
do i2 = nl,n2 

write(6,102) il , i2,Hr2(il , i2) 
enddo 

enddo 

stop 

101 formatC H(',i2,',x) = ',fl8.15) 

102 formatC H(' ,i2, ' , ' ,i2, ' ,x) = ',fl8.15) 
end 

9 Numerical examples 

In Fig. HO, we plot real and imaginary parts for a selection of irreducibe HPLs over the interval [—5; 5]. The 
remaining irreducible HPLs can be obtained by simply interchanging x — > —x, according to Eqs. (|5.1|),(pT2|). 
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10 Summary 

In this paper, we have described the routine hplog, which evaluates the harmonic polylogarithms up to 
weight 4 for arbitrary real arguments. The evaluation is based on a series expansion in terms of appropriately 
transformed expansion parameters for small values of the argument. The evaluation for large arguments is 
based on transformation formulae, relating HPLs of different arguments. The algorithms used and described 
here can be extended to higher weights without further modification, requiring only the values of the HPLs 
in x = 1 to be known. 
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Weight 


HPL 


HPL at x = 1 


w = 1 


H(l;a;) 
H(0;a;) 

H(-l;x) 


oo 



ln2 


w = 2 


H(0,l;x) 

H(0,-l;x) 
H(-l,l;x) 


tt 2 /6 

7T 2 /12 

7r 2 /12-ln 2 2/2 


w = 3 


H(0,0,1 

H(0,1,1 

H(0,0,-1 

H(0,-1,-1 

H(0, -1,1 

H(0, 1,-1 

H(-l,-l,l 

H(-l,l,l 


x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 


C3 

Cs 

3Cs/4 

Cs/8 

-7r 2 ln2/4 + 13C 3 /8 

7r 2 ln2/4 - Cs 

-ln 3 2/6 + C3/8 

-7r 2 ln2/12 + 7C 3 /8 + ln 3 2/6 


w = 4 


H(0, 0,0,1 

H(0, 0,1,1 

H(0, 1,1,1 

H(0, 0,0,-1 

H(0, 0,-1,-1 

H(0, -1,-1,-1 

H(0, 0,-1,1 

H(0, 0,1,-1 

H(0, -1,0,1 

H(0, -1,-1,1 

H(0, -1,1,-1 

H(0, 1,-1,-1 

H(0, -1,1,1 

H(0, 1,-1,1 

H(0, 1,1,-1 

H(-l, -1,-1,1 

H(-l, -1,1,1 

H(-l, 1,1,1 


x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 
x) 


7T 4 /90 

tt 4 /360 

7T 4 /90 

7tt 4 /720 

-7r 2 ln 2 2/12 - tt 4 /48 + 7C 3 ln2/4 + ln 4 2/12 + 2Li 4 (l/2) 

7r 2 ln 2 2/24 + tt 4 /90 - 7C 3 ln2/8 - ln 4 2/24 - Li 4 (l/2) 

-7r 2 ln 2 2/12 - tt 4 /180 + ln 4 2/12 + 2Li 4 (l/2) 

-197r 4 /1440 + 7C 3 ln2/4 

tt 4 /480 

7r 2 ln 2 2/24 - tt 4 /80 + ln 4 2/12 + 2Li 4 (l/2) 

-7r 2 ln 2 2/4 - 7tt 4 /720 + 21C 3 ln2/8 

57r 2 ln 2 2/24 + 7tt 4 /288 - 21C 3 ln2/8 - ln 4 2/12 - 2Li 4 (l/2) 

-11tt 4 /720 + ln 4 2/8 + 3Li 4 (l/2) 

-7r 2 ln 2 2/8 + 7tt 4 /288 - ln 4 2/8 - 3Li 4 (l/2) 

7r 2 ln 2 2/12 - tt 4 /80 + 7C 3 ln2/8 + ln 4 2/24 + Li 4 (l/2) 

7r 2 ln 2 2/24 + tt 4 /90 - 7C 3 ln2/8 - ln 4 2/12 - Li 4 (l/2) 

tt 4 /720 - C3I112/8 + ln 4 2/24 

Li 4 (l/2) 



Table 1: List of irreducible HPLs chosen in the numerical implementation, and their values for x = 1. 
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Figure 1: Plots of irreducible HPLs in the interval [—5; 5]. Irreducible HPLs that follow from the above by 
x — > — x are left out. Solid line: real part, dashed line: imaginary part. 
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